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We study thermodynamics of the 3D Hubbard model at half filling on approach to the Neel 
transition by means of large-scale unbiased Diagrammatic Determinant Monte Carlo simulations. 
We obtain the transition temperature in the strongly correlated regime, as well as temperature 
dependence of energy, entropy, double occupancy, and the nearest-neighbor spin correlation function. 
Our results improve the accuracy of previous unbiased studies and present accurate benchmarks in 
the ongoing effort to realize the antiferromagnetic state of matter with ultracold atoms in optical 
lattices. 



I. INTRODUCTION 

The Hubbard modeP of interacting fermions in a solid 
is a centerpiece of modern condensed matter physics. It is 
conventionally defined by restricting the motion of elec- 
trons in a crystalline solid to a single band, and simpli- 
fying the screened long-range Coulomb interactions be- 
tween electrons to an on-site repulsion: 

(xy)(T X XCT 

(1) 

where a =t, i, Cxo- creates a fermion on a site x, ~ 
c](.g.Cx(T, the summation in the first term runs over the 
nearest-neighbor sites of the simple cubic lattice, t is 
the hopping amplitude and ?7 > is the onsite re- 
pulsion. Though remarkably simple in appearance, the 
model has been used to study a wealth of intriguing 
quantum many-body phenomena that are due to elec- 
tron correlations in solids, such as interaction-driven 
insulator^, quantum magnetism and high-temperature 
superconductivitji^. However, despite more than a half 
of century of intensive investigation, the physics of the 
model is still not completely understood. 

The most challenging yet the most interesting regime 
is the intermediate regime with interaction comparable 
to half the bandwidth U ~ zt, where z is the number of 
nearest neighbors of a site. This regime offers no small 
parameter to start a controllable analytic theory. Fur- 
thermore, exact analytic solutions are only accessible ei- 
ther in one spatial dimensiorP or an infinit^ number of 
spatial dimensions. Substantial progress has been pos- 
sible with the development of efficient quantum Monte 
Carlo (QMC) methods (for a recent review, see Ref. ^ 
accompanied by advances in computer technology. Al- 
though for generic bosonic systems virtually any equilib- 
rium property can nowadays be calculated by QMC with 
a controlled high accuracy^, systematic-error- free simula- 
tions of correlated fermions have been limited to a hand- 
ful of special cases due to the negative sign problerrP. 
The sign problem manifests itself as an exponential scal- 



ing of the simulation time with the system volume and 
inverse temperature making it practically impossible to 
obtain any reliable information about the system in the 
thermodynamic (TD) limit. Although in some cases the 
sign problem can be completely eliminated by choosing 
a system-specific representation, its general solution is 
almost certainly not possibl^. 

A major step toward understanding strongly corre- 
lated systems has been the experimental realization of the 
Hubbard model with ultracold atomic gases loaded into 
optical lattices (for recent reviews see Refs. M and [TUl) . 
These systems offer substantial control over the Hamil- 
tonian. As a result, these experiments can serve as em- 
ulators of quantum many-body systems, which allow the 
accurate study of a given model in a range of param- 
eters inaccessible by analytic and numeric techniq ue J^. 
The recent experimental observation of Mott physicJ^^HSl 
in the Hubbard model is a major milestone along these 
lines. The next crucial step would be a realization of the 
antiferromagnetic (AFM) transition and the Neel state in 
the Hubbard model, which requires a substantial effort 
in reaching lower temperatures, controlling the equilibra- 
tion rates, developing new probing techniques, etc. In 
addition to these inherent challenges, there is also a fun- 
damental problem related to thermometry in ultracold- 
atom systems since they are insulated from the environ- 
ment. Isolated systems such as these are thus, by their 
nature, better characterized by entropy, rather than tem- 
perature. Moreover, probes have to be calibrated in the 
relevant regime and the results obtained with the setup 
validated against available benchmarks. For these pur- 
poses, reliable and accurate numeric results for fermionic 
systems are indispensable because they ultimately al- 
low a full quantitative understanding, as was recently 
demonstrated by the example of a bosonic optical-lattice 
emulatoi'^. 

This work provides reliable benchmarks for the real- 
ization of the Neel state in optical lattices as well as for 
new theoretical methods. We focus on the special case 
of the half-filled {{riy^a) = 1/2, or equivalently /i = U/2) 
three dimensional (3D) Hubbard model (IT]) on a simple 
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cubic lattice. The case of half filling is special due to the 
SU{2) symmetry of the Hamiltonian, which is ultimately 
broken by the Neel state, making magnetism the leading 
instability. In the limit U/t ^ 1, the effects of the in- 
teraction can be studied perturbatively. In the opposite 
limit, U/t > 1, Eq. Q reduces to the AFM Heisen- 
berg model with J ^t^/U. While there is no doubt 
that in the strongly correlated regime the ground state of 
the half-filled model Q is AFM, mapping out the finite- 
temperature phase diagram and studying thermodynam- 
ics of the system is extremely challenginJSHUl, Xo this 
end, we employ the unbiased continuous-time determi- 
nantal diagrammatic Monte Carlo (DDMC) techniqu^l^!^ 
which produces numerically exact (up to a known statis- 
tical error bar) results for a finite-size system, and which 
is free from the fermionic sign problem at half filling on 
bipartite lattices, allowing a reliable extrapolation of re- 
sults to the thermodynamic limit. 

We use DDMC to determine several critical properties 
of the Hubbard model with high control and accuracy 
and compare with high-temperature series expansionPSl 
(HTSE) where possible. We study the range of on-site 
interaction 4 < [//t < 8, where the critical tempera- 
ture Tfq of the AFM transition is expected to reach its 
maximum^- . Results in this regime are vital to op- 
tical lattice emulator efforts offering experiments their 
best chance of observing the AFM phase. We obtain the 
critical temperature Tjv and compute important thermo- 
dynamic properties of the model — the energy and the 
entropy, as well as two optical lattice observables: dou- 
ble occupancy and the nearest-neighbor spin-spin corre- 
lation function — in the paramagnetic phase as a func- 
tion of temperature down to T/y. A number of previous 
works using different unbiased approaches studied T^r (by 
the determinant quantum Monte Carlo (DQMC^f^ and 
the Dynamical Cluster Approximation (DCA^j^ and the 
thermodynamic properties in question (by DDMCP^a 
combination of DCA and DDMCP, and by DQMCPSl) 
in this regime. Our work improves the accuracy of the 
previous results at half filling and provides the most accu- 
rate to date values of the critical temperature T/v and the 
entropy at the critical point Sn (summarized in the Ta- 
ble. |l]). Accurate knowledge of Tjy is particularly impor- 
tant for determining the critical entropy Sn since S{T) 
is a steep function near the transition, so that the error 
bar of Sn is mainly due to the uncertainty in Tjv- The 
values of Sjq are required for an experimental realization 
for the AFM state. As was noted in Ref. [T71 close to the 
transition, the temperature dependence of the nearest- 
neighbor spin-spin correlation function (S'^^x+eJ {^i is 
the unit vector in the direction i) is significantly more 
pronounced than that of the double occupancy (?T.xt"-x^) 
of a lattice site making mea surem ents of the nearest- 
neighbor spin-spin correlationJ^^^^ more suitable for ac- 
curate thermometry in this regime. Our results for the 
spin-spin correlation function can be used for calibration 
of such a thermometer. 

For over a decade, the DQMC study of the phase di- 



agram of the half-filled 3D Hubbard model by Staudt et 
(j^jH! ]^a,s been the main reference for T/y in the corre- 
lated regime. Representing the state of the art at that 
time, Ref. IT provides a comprehensive comparison of 
the DQMC data for the Neel temperature with those of 
preceding QMC simulations and approximate theories, 
e.g., DMFT'^. We shall not reproduce this comparison 
here and refer the reader to a review of results for T^r 
predating Ref. IT^ 

The simulation method of Staudt et al. is based on a 
discrete Hubbard-Stratonovich decoupling in the Hub- 
bard interaction term which requires discretizing the 
imaginary-time interval Q < t < 13 — 1/T into a finite 
number of steps of size Ar thereby introducing a system- 
atic error. Hence, in addition to the standard extrapola- 
tion of the finite-size DQMC data to the TD limit, one 
has to perform an extrapolation with respect to Ar 0. 
Such a double extrapolation is rather laborious. In prac- 
tice Ar is usually fixed at a value which is large enough 
to allow efficient simulations, yet, according to Ref. [TH 
such that "the results are not significantly affected by the 
extrapolation Ar 0" . In the absence of an explicit ex- 
trapolation, the degree of control over systematic errors 
can be questioned. This is where our approach is a ma- 
jor improvement over that of Staudt et al. The DDMC 
technique is formulated directly in continuous imaginary 
time. Therefore finite-size corrections are the only sys- 
tematic error we have to eliminate in our approach. The 
cost of the absence of the additional systematic error is 
the computation complexity of the DDMC, which scales 
as [(3U]^Vt^ versus the linear in [/3/Ar]r2'^ scaling of the 
DQMC of Ref. El for lattices with = sites. As a 
result, we are limited by the values of interaction U < 8, 
whereas Staudt et al. were able to study the model up 
toU = 12. 

We have also been able to identify the transition it- 
self with a considerable improvement in accuracy over 
Ref. [131 Long-range AFM order in the system causes 
divergence of the magnetic structure factor 

where Sz(x) = (rix-f — nx^)/2, and Q = (7r,7r,7r) is the 
AFM wave vector, so that S{Q)/L^ is related to the mag- 
netization m in the TD limit, limi_>.oo S'(Q)/i'^ — m^. 
In Ref. [T31 the transition temperature is found as the 
point at which m starts to noticeably deviate from zero, 
while the magnetization is obtained from a finite-size ex- 
trapolation of S{Q) with respect to L — oo. Here the 
scaling power of the finite-size corrections was used as a 
fitting parameter. An additional (indirect) probe of the 
transition used by Staudt et al. is the peak in the depen- 
dence of the specific heat on temperature. In contrast, we 
use a much more accurate method to determine the crit- 
ical point. We find the critical point using a finite-size 
scaling analysis of S'(Q) in combination with the tech- 
nique of Binder crossings'^^. This approach allows us to 
get a reliable and accurate value of Tjv by making use of 
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the known scahng law of the magnetic structure factor at 
the critical point, S'(Q) oc L^^^, where rj is the anoma- 
lous dimension of the 3D Heisenberg universality class, so 
that the quantity S{Q)L^~^^ becomes scale-invariant at 
the critical temperature up to non-universal corrections, 
which we also take into account. 

In the context of previous calculations of the Neel tem- 
perature, it is worth noting the work by Kent et al^^, 
where Tat was found from calculations using DCA. Al- 
though this work did not lead to any improvement of the 
precision claimed by Staudt et al, it was shown that the 
DCA can be used to determine the critical point with a 
controlled accuracy based on significantly smaller clus- 
ters of only ~ 50 lattice sites versus up to ^ 1000 sites 
in Ref. [TJ] and in this work. In addition, good agree- 
ment with the values of T/v of Ref. [TJ] suggested that 
potential systematic error of the At — >■ extrapolation 
in Ref. TTis likely to be small. In the range of interaction 
U/t > 6, our results for T^r agree within the error bars 
with those of Refs. [HI and \TS\ implicitlv confirming this. 
At smaller U, however, we find somewhat lower values of 
T/v. Moreover, at U/t — 4, being unable to reach signif- 
icantly low temperatures to accurately infer the critical 
point, we can only claim an upper bound from our finite- 
size-scaling analysis, Tff/t < 0.17, which is already lower 
than the values claimed in Refs. "T?' and 15. The reason 
for the discrepancy is likely to be the long-range char- 
acter of correlations at smaller coupling, which can be 
missed in simple finite-size extrapolation schemes based 
on data for insufficiently large systems. 

The thermodynamics of the Hubbard model near the 
Neel transition in connection with its experimental re- 
alization has been the focus of a num ber of recent 
studiePHUmHU. The DMFT resultPl22l emphasize the 
role of double occupancy in detecting the build up of 
AFM correlations. However, its dependence on temper- 
ature near the Neel transition is relatively flat in the 
regime where Tn is maximal, as observed in unbiased (ex- 
trapolated to the TD limit) DCA calculations^^. Fuchs et 
al\^ obtained the energy and the entropy down to Tjv in 
this regime as well as the equation of state away from half 
filling, which allowed to get an estimate of the entropy at 
the transition in a realistic harmonically trapped system. 
A studjff^ of the system using the same DQMC method 
as Ref. [13] arrived at an agreement with the results and 
conclusions of Fuchs et al. . 

The DDMC simulation method used in this work is 
not capable of capturing thermodynamics of the Hubbard 
model away from half filling due to the pronounced nega- 
tive sign problem. However, exactly at half filling, it has 
certain advantages over DCA and DQMC. As discussed 
above, in order to claim unbiased results in the TD limit 
within DQMC one has to resort to a double extrapola- 
tion. At 0, L — >■ cx). DDMC is also formulated directly 
in continuous tim^l^l^ \)ixt the accessible system sizes (e.g., 
10 X 10 X 10 for [/ = 8, T > Tn) are typically larger 
than those of DQMC (8 x 8 x 8 in the same regim^^^ 
allowing a more reliable extrapolation to the TD limit. 



Modern efficient solvers for DCA also work in continuous 
time, but the cluster sizes amenable to simulation in the 
regime of interest are typically less then 100. However, 
in DCA the clusters are embedded in a self-consistently 
defined medium, which largely improves the convergence 
to the TD limit. In practice the finite-size dependence 
for the accessible clusters is notably larger than that in 
DDMC^-^ Hence, at half fiUing and U/t ^ 8, where Tn 
and Sn are expected to reach their maxima, DDMC cur- 
rently allows to obtain the most reliable benchmarks for 
the thermodynamic quantities of interest. 

The paper is organized as follows. In Sec.|llj we discuss 
the simulation method outlining the general formulation 
of the DDMC technique and its application to calculat- 
ing the specific observables in question. 
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Section 

concerned with determining the temperature of the Neel 
transition. Sec. |IV| describes the thermodynamics near 
Tn- Here we discuss the extrapolation of the observables 



to the TD limit (IV A I, the determination of entropy 



(IV C), and thermometry near Tn (IV D). We summa- 



rize the results in Sec. |Vj The Appendix contains tables 
of the obtained numerical data for the entropy, energy, 
double occupancy, and spin-spin correlation functions as 
a function of temperature. 



II. METHOD 

We first rewrite the Hubbard Hamiltonian ([T]) to a form 
suitable for numerical simulations by mapping the repul- 
sive model ([!]) to an attractive model by a particle-hole 
transformation^. We use the fact that the simple cubic 
lattice is bipartite and can be split into two interpene- 
trating sublattices A and B, so that the hopping term 
in ([T]) only connects sites belonging to different sublat- 
tices. Then we introduce the hole operators for the f- 
componcnt: 



Cx-f , xe A 
-Cxt , X e S 



(3) 



This way, Eq. ([I]) becomes at half filling 

iJ = -t ^ a'l^aya + h.c. - mx^TOx; 

(xy)cr X 



fi'^m^^-^n , (4) 



where ax; = Cx;, 



a^g-Cxcr is the number operator 



for the attractive model, fi' = —U/2 as appropriate for 
half filling, and ft = is the total number of sites. 

Since we only consider half filling, (nix^) — (tIx^) — 
1/2, we follow Ref. ISTI and shift the chemical potential 
according to fi' ^ ^' + aU : 



H = Ho + Hi 



a^-^]un 



(5) 
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with 



Hq — -t ^ aitTOycr + h.c. — + all) ^ i 

(xy>Cr X(T 



(6) 



At half fiUing the choice a = 1/2 is optimal because 
it leads to the minimal computational complexity of 
the simulations (see below), and we use this value of a 
throughout. 



A. Diagrammatic Determinantal Monte Carlo 
method 

The DDMC algorithm works with the weak-coupling 
expansion series for the finite-temperature partition func- 
tion for the Hubbard model The latter reads 

oo „ p 

Z^ZoY.U" J2 / []2?(xiTi;...;Xprp) 

p=0 xi...Xp "'0<^i<---<^p</3 j = l 

(7) 



where 

V{xiTi; . . . XpTp) = ( n (™t(xj^i) - a) (miiXjTj) - a)\ 

(8) 

Here (3 is the inverse temperature, Zq ~ TrTexp {—I3Hq) 
is the unperturbed partition function, T denotes time 
ordering, and ((. . . ))o = Tr [T(. . . ) exp (-/37?o)] /Zo de- 
notes the thermodynamic average with respect to the un- 
perturbed Hamiltonian Hq. 

Equations ([7])-(|8]) generate the standard Feynman di- 
agrams: There are (pl)^ diagrams of order p, which can 
be represented graphically as a collection of p vertices 
connected by single-particle propagators. Summing over 
all the interconnections for a fixed vertex configuration 

5p = {(xjrj),j = 1,. (9) 

equation (|8| takes the form 

V{Sp)^\detA{Sp)W (10) 

where A(iSp) are p x p matrices with matrix elements 
given by = l,...,p) 

A,,{Sp) = G(")(x, - Xj, - Tj) - aS,, (11) 

(since we only consider an unpolarized system, we omit 
the spin index for A-s and G^^^-s), and G^") being the 
free-particle Green's functions, 



Since Hq, Eq. ([6]), is diagonal in momentum space, 
Hq = ^[ck - A^' - ctU]al^a]^a, 

kcr 

Sk = -2t J2 cosih), (13) 



1=1,2,3 



the free propagators are calculated by the Fourier trans- 
form 



G(°)(r,r)=^G(")(k,r)e-*'^^ 



G(°)(k,T) = - 



> 0, (14) 



[1 + e-l3{^k-i-'-'-aU)^ 
G(")(k,-r) = -G(0)(k,/3-T), 

and tabulated before the start of the simulation. 

The series ([?]), ([lO])-([Tl]) serves as a basis for a DDMC 
simulation of the Hubbard model ([6]): We set up a ran- 
dom walk in the space of the vertex configurations 5p, Eq. 
([9]), using the Metropolis algorithnPH with the weights 
proportional to I?(5p), Eq. (10 1. Since the technique it- 
self is detailed elsewher G^^ here we only briefly discuss 
the specific details of the present implementation. We 
only stress at this point that since all the terms in the 
series ([7]),([To])-( 11 ) are positive definite we completely 
avoid a sign problem. 

The simplest updating strategy for DDMC simula- 
tions consists of adding (p — > p -f 1) and removing 
(p — !■ ]3 — 1) interaction vertices at random positions in 
x and T to/from a vertex configuration Sp. However, 
at half filling and with a ~ 1/2 the series only contains 
even order terms, hence we employ rank-2 updates, where 
p — >■ p ± 2. Using the Woodbury-type formulas, both 
rank-1 and rank-2 updates can be performed in 0{p^) 
operations. We note that for a 7^ 1/2 both even and odd 
terms are present even at half filling. In this sense the 
choice of a = 1/2 is optimal. 



B. Observables 

The general method for calculating observables in the 
DDMC simulations uses the standard technique of Monte 
Carlo estimators: for an observable O we define an esti- 
mator Q''*^''(5p) such that the average of the latter over 
the vertex configurations generated by the MC process, 
(Q('^)(5p))]yjQ, converges to the thermal average (O), 
where 



(O) = Z-^Tr [Oe-'^"] 



(15) 



Below we explicitly list estimators for useful observables: 
a. Filling fraction The thermal average of the filling 
fraction rUa- for the spin projection a is given by 



G(0)(x. - X,, r, - Tj) - -(ax, (n)4, (r,)>o. (12) 



{ma) = (4^ax 



(16) 
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hence the corresponding estimator 



dot B(5p; xr) 
det A(5p) 



(17) 



where A(iSp) is a p x p matrix (111 and B(5p;xr) is a 
(p + 1) X (p + 1) matrix with an extra row and a column 
corresponding to the extra creation and annihilation op- 
erators in ( |16[ ). Here x and r are random positions in 
space and time, respectively. 

Notice that for a half filled model and for a = 1/2, the 
MC average of the second term in (17) must equal zero. 



We use this fact to check if a simulation has equilibrated. 

b. Kinetic energy Calculating the kinetic energy for 
the model ^ requires evaluating the average of {al.^aya), 
which differs from Eq. ( 16 1 only in that the creation an- 



nihilation is shifted in space with respect to the creation 
operator while in Eq. ( 16 1 both operators reside on the 



same lattice site. The estimator for the kinetic energy 
(Hn) is theiP 



detB(5p;x,y,T) 
det A(5p) 



x2dL-^ 



(18) 



where the matrix B(5p; x, y, r) differs from B(5p; xr) by 
the fact that the creation operator in ( 18 ) is shifted in 



space with respect to the annihilation operator. The ex- 
tra factor of 2 accounts for the summation over a —'[, J, 
and dL^ is the number of bonds of a lattice with sites 
and periodic boundary conditions (PBC). 

c. Interaction energy Since the series ([7| is noth- 
ing but an expansion in powers of Hi , the corresponding 
estimator is readily obtained by a standard trick of con- 
sidering the Hamiltonian Hq + XHi and differentiating 
with respect to A . The result iJ^^ 



P 



(19) 



d. AFM structure factor Calculating the AFM 
structure factor ^ in the particle-hole transformed 
model ([6| requires calculating two independent equal- 
time density-density correlation functions 



(20) 
(21) 



The estimator for the equal spin density correlation 
function g-^^ for x 7^ y is given by 

.^det^,^^^ (22) 
det A(6p) 

where B2(iSp; x, y, r) is a (p -I- 2) x (p -I- 2) matrix with 
two extra rows and columns corresponding to the extra 
creation and annihilation operators in Eq. ( pO| . For x = 
y, g^^{x = y) {ml^] 
model. 



which equals 1/2 for a half- filled 



To build the estimator for Eq. ( 21 1 we proceed similarly 



to (19 1. The resulting estimator for (21) is 



Q(st^)(5p) = 



p detB(5p;y,r) 
pun detA(5p) 



(23) 



where B(5p;y,r) is a p x p matrix constructed by se- 
lecting a random vertex, (xr), from a configuration Sp 
and moving the corresponding row of the matrix A(5p) 
to (yr) while leaving the corresponding column at (xr) . 



III. CRITICAL TEMPERATURE 

In the paramagnetic phase, T > T/v, S{Q) scales to 
zero exponentially as L — ?> 00. In the AFM phase, on the 
other hand, S{Q) — > as i ^ 00, where m is the sub- 
lattice magnetization. Right at the critical temperature 
5(Q) cx 1/L^+'', where 77 is the anomalous dimension. 

In order to locate the transition temperature we thus 
use the standard finite-size scaling (ESS) ansatd^^ 



S{Q)L 



i+i) _ 



f{L/0{l + cL-^ + ■■■) , (24) 



where ^ is the correlation length which diverges at the 
transition as ^ cx jT — Tn\~'^ , f{x) is a real-valued func- 
tion which tends to a finite constant as a; — > 0, and the 
corrections in brackets arise from the leading irrelevant 
operators (dots represent the higher-order corrections). 
Here the exponent w is universal, but the amplitude c is 
not. Assuming criticality in the 3D Heisenberg univer- 
sality class, we take 77 « 0.037, v ~ 0.71 and w « 0.^^ 



The basic idea of using Eq. (24) for the FSS is as 
follows: if the corrections-to-scaling (the 2nd term in 
brackets in Eq. ( p4| ) were not present, 5(Q)i^+'' would 
be scale independent at the transition point, so that 
performing the simulations at a series of system sizes 
Li > L2 > ... and plotting 5'(Q)L^+'' versus temper- 
ature, one would observe that all the curves intersect at 
the same point, T — Tpf. This is what we indeed ob- 
serve for (relatively) large values of U/t: for U > 6t our 
MC results are consistent with c = in Eq. (24) within 
statistical errors, see Fig. [l] 

We find that the corrections-to-scaling become more 
pronounced with decreasing U jt: For U = 5t we find a 
clear evidence of the shift of the pairwise crossings to- 
wards lower temperatures, see Fig. [2] Since simulating 
larger system sizes is not an option, we employ Eq. (24) 
including corrections to scaling. The most straightfor- 
ward way is to follow the evolution of the pairwise cross- 
ings with the system size. Indeed, expanding Eq. (24) 
around the crossing of S{Q)L^~^^ at system sizes L = Li 
and L = L2 up to the terms linear in T — T/v we find (cf 
Ref.lini): 



TliM -Tn ^ const X g(Li, L2) 



where 



giLi,L2 



{L2/L1 



1 - (Li/L2)i/^ 



(25) 



(26) 
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FIG. 1. (Color online.) Finite size scaling for Tjv at (7 = 8t 
by Binder crossings analysis. Points are Monte Carlo results, 
lines are linear fits. The uncertainty for the Neel temperature 
is estimated conservatively by varying the Monte Carlo points 
within their respective errorbars. 



FIG. 3. (Color online.) Scaling of corrections to the Binder 
crossings, Eq. (25l-(26l at [/ = 5t. See text for discussion. 



We perform the linear fit of the series of crossings T^^ 
versus g{Li,L2). Then the intercept of the best-fit line 
yields the Neel temperature. This procedure is illustrated 
in Fig.[2l 
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FIG. 2. (Color online.) Finite size scaling for Tjv at [/ = 5t 
by Binder crossings analysis. Points with errorbars are Monte 
Carlo results, solid lines are linear fits. FSS procedures based 
on Eqs. (25l-(26l and ( |27[ l result in the estimates Tjv/t = 
0.2175(44)and Tjv/t = 0.2211(26), respectively. See also Fig. 
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An equivalent procedure has been suggested in Ref. [35l 
Again, one expands Eq. (24 1 up to the linear order in 
T — T/v, which leads to 



ai{T-TN)L^/'') (l + cL 



which is then used as a four-parameter ansatz for a single 
nonlinear fit. A priori, fitting procedures based on (25) 
and (27) are equivalent and indeed produce consistent re- 
sults, as illustrated in Fig. [2j We stress at this point that 
using (27) requires judicious choice of the temperature 



range for fitting: including Monte Carlo points at too 
high temperatures and/or too small system sizes tends 
to significantly skew the fit results. In the following we 
therefore quote the T^r-s obtained using Eq. ( 25 ) . 



For U = 4t, we find the corrections-to-scaling to be 
larger than those for U = 5t, see Fig.|4] In fact, with the 
accessible systems sizes we are only able to put an upper 
limit on the Neel temperature, T/v < O.lTt. From Fig. |4] 
it is clear that L = 6 and possibly even L = 8 are simply 
too small and need to be discarded from the finite-size 
scaling analysis. 



Our results for the dependence of the Neel tempera- 
ture on U are summarized in Table |T] and Fig. [5] It is 
instructive to compare our estimates to the previous un- 
biased calculations from the literature. While for U/t = 6 
and 8 our estimates agree with and are more accurate 
than previous estimates from QMCJ^ and DCA^. For 
smaller values of U/t our estimates are systematically 
lower. The discrepancy can be traced back to the FSS 



procedure which includes corrections-to-scaling, Eq. (24 1 



(27) 



if we were to discard the corrections and identified the 
Binder crossings oi L — 6 and L = 8 with the Neel tem- 
perature, such estimates would have agreed with Refs. [14] 
and US] We therefore conclude that the estimates of Tjy 
presented here are more accurate than results reported 
to date. 
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IV. THERMODYNAMICS 



5(Q)L 




0.14 016 0J8 020 022 024 f/t 



FIG. 4. (Color online.) Finite size scaling for Tjv at U = it by 
Binder crossings analysis. In this case we not able to reliably 
extract the Neel temperature and can only provide an upper 
limit, Tn < O.nt. Notice that the crossing of L = 6 and 
L = 8 is clearly outside of the range of applicability of either 
(251 or (pTl. 



U/t 




Sn 


4 


< 0.17 


< 0.17 


5 


0.2175(44) 


0.135(25) 


6 


0.300(5) 


0.305(35) 


8 


0.3325(65) 


0.33(3) 



TABLE I. Neel temperatures and entropies. See text for dis- 



cussion 
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FIG. 5. (Color online.) Comparison of estimates for Tjv 
by different unbiased approaches. Also shown is the strong- 
coupling limiting behavior, Tjv = 3.83t^/f/. See text for dis- 
cussion. 



A. Extrapolation to the thermodynamic limit 

The dependence of local observables on the size L of 
the system with PBC is complicated'^^ by oscillations be- 
tween the results for even and odd values of L/2. The 
issue is illustrated in Fig. [6j where the energy per particle 
of the half-filled non-interacting system (J7 = 0, /.t = 0) 
with PBC is plotted versus L^^ up to a large system 
size (L = 52) for different temperatures. The TD-limit 
value is approached from above by the data for even 
i/2 and from below by those with L/2 odd. These are 
the well-known "shell" oscillations'^^ caused by whether 
or not the spectrum of the finite system has states v 
with the energy E^, within a range much less than T 
from the Fermi level, \Ey — /x| <C T. In the exam- 
ple of Fig. |6] the states are classified by the momenta 
fc — {ki,k2,k^), hi — 2Trni/L with integers taking 
the values n, = . . . , -1, 0, 1, . . . , L/2 - 1. When 

{rii] ~ L/A, which is only possible if L/2 is even, the 
state k = (7r/2, 7r/2, 7r/2) is exactly at the Fermi level; 
it's occupation is 1/2 ("open shell"), but it gives no con- 
tribution to the average energy. Hence, if L is not large 
enough so that the spacing between the levels is larger 
than T, the average energy per particle of the system with 
a closed shell {L/2 odd) is systematically lower than that 
of the system with an open shell {L/2 even) due to the 
difference of the number of states below the Fermi level, 
L^ /2 and L^/2 — 1 correspondingly. However, for any 
given temperature T there is a system size — (T) 
such that for L > L^, the number of states with the en- 
ergies \E,j — fi\ < T becomes large removing the distinc- 
tion between even and odd L/2. In the free-particle case 
of Fig. [6j the convergence to the TD limit dX L > L^, 
is extremely fast (exponential) with L*(r — 0.5) ^ 10, 
L.,{T = 0.3) - 16, and i*(T = 0.1) - 40. 

At the finite values of U studied here, the param- 
agnetic phase should be described by a Fermi liquid 
in the limit of T <C Ep- In this regime, the total 
energy is a functional of occupation numbers of non- 
interacting quasiparticles. Therefore, the system-size de- 
pendence of energy is expected to be proportional to that 
of the non-interacting system, at least for large enough 
systemPSI. This implies a TD-limit extrapolation in the 
ioTiaE{L) = E{oo) + C[Eo{L)- Eo{oo)]+ g{L) suggested 
in Ref. |3H1 where C is a constant, Eq{L) is the energy of 
the corresponding non- interacting system of size L, and 
g{L) is an unknown in our case function. One can ex- 
pect that \g{L)\ < \C[Eo{L) - £^o(oo)]| for sufficiently 
large L. Given that our simulations are limited to sys- 
tem sizes of up to L ^ 10, the validity of this condition 
is not guaranteed a priori. An example of such an ex- 
trapolation with g{L) — for a typical set of parameters 
([/ = 8, T = 0.3875) is shown in Fig. [t) which suggests 
that the additional corrections given by g{L) should be 
small for L > 6. Nonetheless, we claim the TD-limit 
value E{L — oo) and its error bar AE{L — )■ oo) con- 




FIG. 6. (Color online.) Dependence of the energy per particle 
of non-interacting fermions at half- filling (described by Eq. ([T]) 
with U — 0, fi = ) on the inverse of the linear system size L. 



servatively as the span between the values at the two 
largest accessible system sizes including their statistical 
error bars (depicted by the horizontal band in Fig.[7]): 



E[L oo] 
minfiJfLma, 



-AE[L 



-2])- 



max{E[Ln,^^]+AE[L^^^],E[L,^^^^2]+AE[L„ 



/2: 



FIG. 7. (Color online.) Example of the dependence of en- 
ergy on the inverse of the linear system size L obtained with 
the periodic boundary conditions (PBC, circles) and using 
the twist-averaged boundary conditions (TABC, triangles) for 
U = 8, T = 0.3875. The solid line is an extrapolation us- 
ing the formula E{L) = E'{cx) + C[Eo{L) - Eo{oo)] (see 
text) with E'{oo) = -0.5965, and C = 0.6. The claimed 
thermodynamic-limit result E{oo) — —0.5960(16) is shown 
by the horizontal band. 



PBC. Then an observable A{L) is obtained by means of 
the integration 



(28) 



and 



AE[L ^ ^] « |£;[L,nax] - E[L„,,^ - 2]|/2 
+ AE[L^,^] + AE[L 

max 



(29) 



and similarly for other local observables. 

As a consistency check for our TD-limit results, we em- 
ploy two other simulation setups, which exhibit different 
system-size dependences, which we detail below. 



1. Twist- averaged boundary conditions. 

Averaging over twisted boundary conditions was found 
in Refs. \S7\ and [5^ to produce exact results for the non- 
interacting system in the grand canonical ensemble and 
to substantially suppress the system-size dependence for 
interacting systems. In this approach one introduces a 
finite phase that particles acquire when they wrap around 
the periodic boundaries. 



Le, 



r2, ■ 



|ri,r2,...), i = 1,2,3, (30) 



where is the unit vector in the direction i and —n< 
Qi < tt; Eq. ( 30 1 with 9i = corresponds to the standard 



AiL) 



1 



Q 



Q 



A&{L)d&, 



(31) 



where A@{L) is a result of the simulation with a fixed 
value of = {©i} and system size L. Thereby, the 
possible momentum values are forced to span the whole 
Brillouin zone. The non-interacting propagators G^^^ sat- 
isfying the condition ( |30[) a re obtained by substituting 
k — > k+0/L in Eq. (|14[). In practice, we perform 
numerical integration on a mesh of 64 points, esti- 
mating the systematic error of integration to be smaller 
than the statistical error of A{L) coming from sampling 
each A@ [L) by Monte Carlo. The results of the calcu- 
lation of energy with the twist-averaged boundary con- 
ditions (TABC) are compared to those for the PBC in 
Fig. [7j The twist-averaging is seen to somewhat reduce 
the finite-size corrections so that the data for L = 8 and 
L = 10 appear to be converged within their error bars to 
the value consistent with the claimed TD-limit result. 



2. L — > oo free propagators. 



In the second approach, we replace the free-particle 
propagators G^"-* = c'f^^ in the diagrammatic expansion. 



G^""* G)x! , thereby completely eliminating the oscilla- 
tions coming from the discreteness of the spectrum at the 



Eq. (11), by those corresponding to the limit L — cx). 



-.(0) 
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l/L 
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dE at fixed volume by the integral 



FIG. 8. (Color online.) Example of the dependence of energy 
on the inverse of the hnear system size L obtained with the 
periodic boundary conditions (PBC, circles) compared to the 
result of a simulation based on free-particle propagators of an 
infinite system (G 



(0) 



G^S, squares) for 17 = 8, T = 0.3875. 



The dashed line is a linear fit yielding E" {oo) = -0.5962(9) 
in perfect agreement with the claimed conservative estimate 
E(oo) = —0.5960(16) shown by the horizontal band. 



expense of giving up the PBC. The only source of system- 
atic error in this case is the finiteness of the volume — still 
given by — confining the distribution of the interaction 
vertices in Eq. (|9|. In this case, the finite-size corrections 
are substantially larger than those of simulations with 
PBC. However, the scaling of these corrections is linear 
in l/L for all the local observables in question, which 
allows to perform a systematic TD-limit extrapolation. 
As an example, such an extrapolation for energy in com- 
parison with the data for PBC is shown in Fig. |8] The 
TD-limit value obtained thereby is in perfect agreement 
with the result of simulations with the PBC. This consti- 
tutes an independent verification of the accuracy of the 
claimed results. 



B. Energy 

Here we present simulation results for the total energy 
per particle for a range of the interaction U in the corre- 
lated regime near the Neel transition. The temperature 
dependence of the energy per particle for U — 8, 6, 5, 4 
is plotted in Fig. [o] along with the results of the HTSE^o 
for orders 2, 4, 6, 8, 10. As seen from the plot, the HTSE 
starts diverging well above the transition point. 



C. Entropy 

The entropy per particle S{T) at a given temperature 
T is obtained from the thermodynamic relation TdS = 



S{T) = S{T,) 



E{T) E{T, 



T 



E{T') 



dT', (32) 



where T* is some temperature at which the entropy is 
known. We choose to be the lowest temperature at 
which the HTSE for the energy obviously converges to the 
TD-limit value E{T^) from the simulation. From Fig. |9j 
we find = 1.8, 2.4, 2.6, 2.6 at C/ = 8, 6, 5, 4 respec- 
tively. Then, the accurate value of S{T^) in Eq. (|32| is 
given by the HTSE, while the integral is done over the 
simulation data after taking the TD limit. Since the de- 
pendence E{T) is slow, we represent it by a piecewise 
linear function and take the integral analytically. The 
systematic error of integration is included in the error 
bars for S{T), but is negligible compared to the error 
propagated from the values of E{T). 

The resulting curves of S{T) for U — 8,6,5,4 are 
shown in Fig. [T0| From these data and our calculation of 
T/v discussed in Sec. |III[ we find the values of the critical 
entropy Sn — S{Tj\j) in the range of U and summarize 
the results in Table|l] The error bars of Sn are dominated 
by the relatively small error of Tjv due to the large slope 
of Sn{T) near the transition. In Fig. 11 we plot lines of 
constant entropy in the (T, U) plane. The latter demon- 
strate that an adiabatic increase of the coupling U can 
lead to either a rise (at S > 0.7) or a fall (at S < 0.7) of 
temperature, although the net effect of the Pomeranchuk 
cooling near Sn is rather small. 



D. Thermometry 

As was shown in Ref. [T7] by means of DCA calcu- 
lations, the nearest-neighbor spin-spin correlation func- 
tion defined as {S^S^_^_^.), which is accessible in present- 
day ultracold-atom experiments, can serve as a sensitive 
thermometer at temperatures near Tn- In contrast, an- 
other routinely measured correlator, the double occu- 
pancy (nx^nx^) of a lattice site, is nearly flat in this 
temperature range making it a rather poor candidate 
for thermometry. This is hardly surprising since the lat- 
ter is concerned with correlations in the charge channel, 
whereas the relevant physics at these temperatures is that 
of developing short-range spin correlations. 

we present our results for (<S'xS'x+e ) ^^'^ 

5,6,5,4. 
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In Fig. 

('^xt^xi) extrapolated to the TD limit at U — 
The obtained values agree within the errors with the 
TD-limit-extrapolated data from the DCA simulations, 
Ref. I17L but our error bars are notably smaller. Our data 
can be directly used for thermometry calibration and de- 
tection of the Neel transition. Note that in the range 
of temperatures Tn < T < 2Tn (the position of Tn is 
depicted by a vertical line with the width correspond- 
ing to the error bar) the spin-spin correlations between 
nearest-neighbor cites rise by a factor of two, with a sub- 
stantial increase of the slope close to Tn- In the same 
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FIG. 9. (Color online.) Energy versus temperature at 17 = 8,6,5,4. The lines represent the results of the high-temperature 
series expansion series of orders 2, 4, 6, 8, 10 labeled correspondingly. 



temperature range the double occupancy varies by less 
than 10%. 



V. CONCLUSIONS 

We presented unbiased results for the 3D Hubbard 
model at half filling near the Neel transition in the 
strongly correlated regime of 4 < [/ < 8, where T^r 
reaches its maximum. We focused on the properties 
of the model near the transition accurately determining 
Tjv and studying the energy, entropy, double occupancy 
(rix^JT-xi) and the nearest-neighbor spin-spin correlator 
(S'^S'x+e ) as functions of temperature and interaction. 
Accurate quantitative understanding of the model in this 
regime is of growing importance in view of the ongoing 
experimental effort to emulate the Hubbard model with 
ultracold atoms in optical lattices, which could ultimately 
allow to study regions of the phase diagram inaccessi- 
ble by unbiased theoretical methods. In particular, this 
could lead to answers of fundamental questions regard- 
ing the nature of superfluidity at finite doping and its 



connection to high-temperature superconductors^^. The 
realization of the Neel state would be a necessary step 
on the way to accessing the region of the phase diagram 
where quantum fluctuations play an important role. Our 
simulations provide the most accurate and controlled es- 
timates of entropy at the critical point to date. These 
entropies, summarized in Table |lj have to be achieved 
in the middle of the trapped cold-atom system to real- 
ize the AFM state. For independent in situ thermom- 
etry in this regime, one can employ measurements of 
the nearest-neighbor spin correlations, which expectedly 
have pronounced temperature dependence near TjVi and 
which can nowadays be addressed either by the use of 
superlattices--' or by lattice modulatiorp2-. In agreement 
with Ref. Ull we did not find the double occupancy to 
display notable temperature dependence in the regime of 
interest. More generally, our results for thermodynamics 
at half filling quantitatively agree with the extrapolated 
DC A data of Ref.ilTiand Ref. il8l with the combined error- 
bars, although the energy and entropy in DCA^*" appear 
to be systematically above our values as well as those of 
DQMC on approach to the critical point. As a result and 
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FIG. 10. (Color online.) Entropy versus temperature a,t U — 8, 6, 5, 4. In each panel, the vertical line shows the position of the 
critical temperature Tjv with its width given by the error bar, while the horizontal dashed lines represent the corresponding 
bounds of the critical entropy Sn listed in Table |l] The rest of the lines represent the results of the high-temperature series 
expansion series of orders 2, 4, 6, 8, 10 labeled correspondingly. 
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FIG. 11. (Color online.) Lines of constant entropy in the T 
vs U plane. 



due to the improved estimate of T/v, our value of Sn at 
U = S (0.33(3)) is below that claimed in Ref. [171(0.42(2)) 
suggesting agreement at the level of two combined stan- 
dard deviations. 

The need for a more precise knowledge of T/v comes 
from the steep temperature dependence of the entropy 
close to the transition. Our results for Tjv improve on 
the earlier studies of Staudt et al}'^, although remain in 
perfect agreement with the latter within the error bars 
everywhere but at U — A, where we were able to find only 
the upper bound for T/v, which is somewhat lower than 
the result of Ref. O The results of our simulations can 
be used as benchmarks for tuning approximate methods 
as well as in developing new unbiased techniques. 

The simulations were carried out on the Brutus clus- 
ter at ETH Zurich. E.K. acknowledges financial sup- 
port of the Fellowship for Advanced Researchers by 
the Swiss National Science Foundation. E. B. grate- 
fully acknowledges the hospitality of Laboratoire de 
Physique Theorique et Modeles Statistiques, where a 
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FIG. 12. (Color online.) Double occupancy (rixt^x^) (squares) and the nearest-neighbor spin correlation function KSxS'x+ej)! 
(circles) versus T at U = 8,6, 5, 4. In each panel, the vertical line shows the position of the critical temperature for a given 
value of U listed in Table |l] with its width corresponding to the error bar. 
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VI. APPENDIX TABLE II. [7 = 8: Energy, entropy, double occupancy, and 

the nearest-neighbor spin-spin correlator as functions of tem- 
perature extrapolated to the TD limit. 



T/t 


E/t 


S 






1.8182 


-0.1434(24) 


- 


- 


- 


1.7 


-0.1855(15) 


0.9408(16) 


- 


- 


1.5385 


-0.2397(26) 


0.9073(22) 


- 


- 


1.45 


-0.2705(28) 


0.8867(23) 


- 


- 


1.4 


-0.2870(31) 


0.8751(26) 


- 


- 


1.3333 


-0.3081(17) 


0.8596(19) 


- 


- 


1.25 


-0.3314(38) 


0.8416(33) 


- 


- 


1.1765 


-0.3515(45) 


0.8250(41) 


- 


- 


1.1 


-0.3802(14) 


0.7998(19) 


0.0765(1) 


-0.0146(4) 


1.05 


-0.3894(35) 


0.7913(36) 


- 


- 


1.0 


-0.4048(28) 


0.7762(31) 


- 


- 


0.95 


-0.4198(22) 


0.7608(27) 


0.0746(3) 


-0.0164(5) 


0.9 


-0.4314(30) 


0.7483(36) 


- 


- 


0.85 


-0.4492(44) 


0.7279(54) 


0.0732(3) 


-0.0193(11) 


0.8 


-0.4607(20) 


0.7140(30) 


0.0727(3) 


-0.0200(3) 


0.75 


-0.4758(19) 


0.6945(29) 


0.0725(2) 


-0.0218(4) 


0.65 


-0.5038(31) 


0.6544(51) 


0.0716(4) 


-0.0250(10) 


0.6 


-0.5157(21) 


0.6354(42) 


0.0717(1) 


-0.0271(2) 


0.55 


-0.5323(11) 


0.6065(26) 


0.0718(3) 


-0.0292(4) 


0.5 


-0.5446(46) 


0.5830(93) 


0.0719(4) 


-0.0311(9) 


0.45 


-0.5647(29) 


0.5407(76) 


0.0714(5) 


-0.0347(10) 


0.3875 


-0.5960(16) 


0.4657(57) 


0.0716(2) 


-0.0410(4) 


0.375 


-0.5967(76) 


0.464(20) 


0.0714(7) 


-0.0424(9) 


0.35 


-0.6240(19) 


0.3887(73) 


0.0719(4) 


-0.0469(7) 


0.335 


-0.6414(19) 


0.3379(59) 


0.0714(2) 


-0.0518(6) 


0.328 


-0.6504(19) 


0.3106(60) 


0.0715(2) 


-0.0537(3) 


0.325 


-0.665(11) 


0.266(36) 


0.0710(7) 


-0.0528(14) 


0.31 


-0.6664(74) 


0.261(24) 
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T/t 


E/t 


S 


(n-xt^x;) 




2.4 


-0.0661(1) 


_ 


0.13908(2) 


-0.00487(3) 


2.2 


-0.1346(1) 


1.0847(1) 


0.13452(3) 


-0.00558(3) 


2.0 


-0.2057(2) 


1.0508(1) 


0.12990(4) 


-0.00639(6) 


1.8 


-0.2787(2) 


1.0124(1) 


0.12533(4) 


-0.00746(5) 


1.6 


-0.3533(2) 


0.9685(1) 


0.12093(4) 


-0.00880(6) 


1.4 


-0.4275(3) 


0.9189(2) 


0.11705(4) 


-0.01041(8) 


1.2 


-0.5003(4) 


0.8628(3) 


0.11380(5) 


-0.01221(13) 


1.0 


-0.5717(6) 


0.7977(6) 


0.11170(7) 


-0.01489(19) 


0.8 


-0.6406(6) 


0.7208(8) 


0.11108(7) 


-0.01833(30) 


0.6 


-0.7054(15) 


0.6276(26) 


0.11298(25) 


-0.02205(40) 


0.5 


-0.7423(15) 


0.5604(35) 


0.11403(29) 


-0.02580(51) 


0.4 


-0.7759(10) 


0.4854(32) 


0.11620(10) 


-0.02885(17) 


0.35 


-0.7962(16) 


0.4312(48) 


0.11685(25) 


-0.03149(26) 


0.325 


-0.8091(17) 


0.3929(54) 


0.11662(40) 


-0.03378(49) 


0.3077 


-0.8251(20) 


0.3423(66) 


0.11521(50) 


-0.03691(87) 


0.3030 


-0.8288(20) 


0.3303(67) 


0.11510(40) 


-0.03800(60) 


0.2963 


-0.8414(50) 


0.288(17) 


0.1134(10) 


-0.04046(34) 



TABLE III. U = Q: Energy, entropy, double occupancy, and 
the nearest-neighbor spin-spin correlator as functions of tem- 
perature extrapolated to the TD limit. 
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T/t 


E/t 


S 






2.6 


-0.0987(6) 


— 


_ 




2.5 


-0.1305(4) 


1.1483(3) 


— 


— 


2.4 


-0.1623(5) 


1.1354(3) 


_ 


— 


2.3 


-0.1957(5) 


1.1212(3) 


_ 


— 


2.2 


-0.2297(7) 


1.1060(4) 


_ 


— 


2.1 


-0.2647(10) 


1.0898(5) 








2.0 


-0.3011(6) 


1.0720(3) 


— 


_ 


1.9 


-0.3371(13) 


1.0535(7) 


_ 


_ 


1.8 


-0.3749(15) 


1.0331(8) 


_ 


_ 


1.7 


-0.4132(16) 


1.0112(9) 


_ 


_ 


1.6 


-0.4517(17) 


0.9879(11) 


_ 


_ 


1.5 


-0.4897(17) 


0.9633(12) 


_ 


_ 


1.4 


-0.5301(9) 


0.9354(7) 


_ 


_ 


1.3 


-0.5687(22) 


0.9069(16) 


_ 


_ 


1.25 


-0.5889(22) 


0.8910(18) 


_ 


_ 


1.2 


-0.6096(14) 


0.8741(12) 


0.1343(1) 


-0.0115(7) 


1.1 


-0.6478(32) 


0.8409(29) 


0.1336(3) 


-0.0121(12) 


1.0 


-0.6886(4) 


0.8020(12) 


0.1328(1) 


-0.0140(1) 


0.9 


-0.7250(36) 


0.7637(40) 


0.1330(5) 


-0.0147(10) 


0.8 


-0.7637(22) 


0.7181(33) 


0.1334(4) 


-0.0161(6) 


0.7 


-0.7996(41) 


0.6701(61) 


0.1348(6) 


-0.0169(15) 


0.6 


-0.8394(15) 


0.6087(45) 


0.1360(3) 


-0.0195(6) 


0.5 


-0.8726(11) 


0.5482(30) 


0.1389(3) 


-0.0203(7) 


0.4 


-0.9110(31) 


0.4626(80) 


0.1411(8) 


-0.0231(12) 


0.35 


-0.9296(17) 


0.4128(81) 


0.1421(4) 


-0.0253(1) 


0.33 


-0.9379(45) 


0.389(14) 


0.1423(12) 


-0.0258(16) 


0.3125 


-0.9404(50) 


0.381(17) 


0.1423(6) 


-0.0266(5) 


0.28 


-0.9535(30) 


0.336(13) 


0.1438(8) 


-0.0262(14) 


0.26 


-0.9665(21) 


0.2884(94) 


0.1423(4) 


-0.0296(6) 


0.25 


-0.9771(19) 


0.2468(78) 


0.1407(10) 


-0.0315(3) 


0.235 


-0.9900(18) 


0.1936(77) 


0.1381(10) 


-0.0361(7) 


0.22 


-1.0016(16) 


0.1424(76) 


0.1362(11) 


-0.0390(7) 


0.21 


-1.0079(10) 


0.1133(53) 







TABLE IV. (7 = 5: Energy, entropy, double occupancy, and 
the nearest-neighbor spin-spin correlator as functions of tem- 
perature extrapolated to the TD limit. 



T/t 


E/t 


S 






2.6 


-0.2218(2) 


_ 


_ 




2.4 


-0.2837(3) 


1.1534(1) 


_ 




2.3 


-0.3164(5) 


1.1395(2) 


_ 




2.2 


-0.3504(3) 


1.1244(1) 


_ 




2.1 


-0.3854(6) 


1.1081(2) 


_ 




2.0 


-0.4213(3) 


1.0906(2) 


_ 




1.9 


-0.4583(4) 


1.0716(2) 


— 




1.8 


-0.4967(4) 


1.0508(2) 


_ 




1.7 


-0.5363(6) 


1.0282(4) 


— 




1.6 


-0.5759(4) 


1.0042(3) 


_ 




1.5 


-0.6174(7) 


0.9774(5) 


_ 




1.4 


-0.6594(5) 


0.9484(4) 


_ 




1.3 


-0.7021(8) 


0.9168(6) 


_ 




1.2 


-0.7445(4) 


0.8829(4) 






1.1 


-0.7875(12) 


0.8455(11) 


0.1553(2) 


-0.0121(6) 


1.0 


-0.8296(17) 


0.8053(17) 


0.1554(2) 


-0.0127(8) 


0.9 


-0.8727(20) 


0.7599(24) 


0.1557(5) 


-0.0138(7) 


0.8 


-0.9152(15) 


0.7098(21) 


0.1562(3) 


-0.0149(5) 


0.7 


-0.9568(29) 


0.6543(42) 


0.1574(6) 


-0.0157(11) 


0.6 


-0.9967(23) 


0.5927(46) 


0.1597(7) 


-0.0164(11) 


0.5 


-1.0390(20) 


0.5155(49) 


0.1613(7) 


-0.0193(8) 


0.4 


-1.0725(44) 


0.441(12) 


0.1654(13) 


-0.0186(18) 


0.3 


-1.1101(19) 


0.333(15) 


0.1665(5) 


-0.0218(7) 


0.27 





_ 


0.1671(9) 


-0.0226(7) 


0.25 


-1.1243(25) 


0.281(12) 


0.1680(7) 


-0.0221(8) 


0.2 


-1.1390(9) 


0.215(11) 


0.1678(4) 


-0.0237(8) 


0.1818 


-1.1462(8) 


0.1775(56) 


0.1657(5) 


-0.0265(5) 


0.16 


-1.1515(31) 


0.147(20) 


0.1639(30) 


-0.0278(34) 



TABLE Y. U = A: Energy, entropy, double occupancy, and 
the nearest-neighbor spin-spin correlator as functions of tem- 
perature extrapolated to the TD limit. 



